// zero sized data


#include <iostream>
	using std::cerr;
	using std::ostream;
#include <limits>
	using std::numeric_limits;
#include <iterator>
	using std::back_inserter;
#include <string>
	using std::string;
#include <algorithm>
	using std::max;
	using std::min;
#include <vector>
	using std::vector;
#include <cassert>
#include <unordered_map>
	using std::unordered_map;


#include <boost/range/algorithm.hpp>
	using boost::range::lower_bound;
#include <boost/range/adaptor/map.hpp>
	using boost::adaptors::map_keys;
#include <boost/range/algorithm/copy.hpp>
	using boost::range::copy;
#include <boost/filesystem.hpp>
	namespace fs = boost::filesystem;
#include <boost/filesystem/fstream.hpp>
	//using fs::ifstream

#include <boost/spirit/include/karma.hpp>
#include <boost/spirit/include/qi.hpp>
	namespace qi = boost::spirit::qi;
//#include <boost/spirit/include/support_multi_pass.hpp>
#include <boost/fusion/include/vector_tie.hpp>                                                                                                 
#include <boost/fusion/include/adapt_struct.hpp>
#include <boost/fusion/include/io.hpp>

#include <boost/format.hpp>
	using boost::format;
//#include <boost/assign.hpp>
//#include <boost/assign/std/vector.hpp> 
#include <boost/assign/list_of.hpp> 
	using namespace boost::assign; // bring 'list_of()' into scope

#include "intervals_per_contig.h"
#include "sort_by_embedded_numbers.h"
#include "sort_by_canonical_chromosomes.h"

#define NDEBUG 1
#ifndef	NDEBUG
#define DEBUG_CERR(text) {std::cerr << text;}
#else
#define DEBUG_CERR(text)
#endif

struct t_parsed_genomic_interval
{
    string  contig;
    int     beg_pos;
	int		end_pos;
};


BOOST_FUSION_ADAPT_STRUCT(
t_parsed_genomic_interval,
(std::string, contig)
(int, beg_pos)
(int, end_pos)
)

// N.B. uses open-close coordinates

namespace intervals_per_contig
{


//!_______________________________________________________________________________________
//!
//! 	t_interval::str
//! 
//! 			for debugging purposes mainly
//!     
//!     \author lg (1/7/2011)
//!     
//!     \return string 
//!_______________________________________________________________________________________
string t_interval::str() const
{
	if (!id)
		return (format("[%1$d-%2$d]") % beg  % end).str();
	else
		return (format("[%1$d-%2$d : %3$d]") % beg  % end % id).str();
}


//!_______________________________________________________________________________________
//!
//! 	out_of_bounds_intervals
//! 
//! 		special magic intervals to indicate out of bounds
//!     
//!     \author lg (1/7/2011)
//!_______________________________________________________________________________________
t_vec_intervals t_interval_search_cache::out_of_bounds_intervals = list_of(t_interval(0,0))(t_interval(numeric_limits<long>::max(),numeric_limits<long>::max()));



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     t_interval_search_cache definitions
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888

//!_______________________________________________________________________________________
//!
//! 	t_interval_search_cache::str
//! 
//! 			for debugging purposes mainly
//!     
//!     \author lg (1/7/2011)
//!     
//!     \return string 
//!_______________________________________________________________________________________
string t_interval_search_cache::str() const
{
	if (pintervals == 0)
		return ("<search cache = 0>");
	else
		return (format("< %1$s - %2$s >") % i_prev->str()  % i_next->str()).str();
}


//!_______________________________________________________________________________________
//!
//! 	t_interval_search_cache::operator++()
//! 
//!     	Increments alternately to flank a interval, or to two intervals flanking a gap
//!     
//!     \author lg (1/7/2011)
//!     
//!_______________________________________________________________________________________
void t_interval_search_cache::operator++()
{
	assert(pintervals);
    // can't increment any more
    if (next_out_of_bounds())
    {
        DEBUG_CERR("    increment_flanking_interval_pair: Already beyond any intervals\n")
        return;
    }

    // 
    // Both point to the same interval,
    // Make them flanking by increment the next iterator
    // 
    if (single())
    {
		++i_next;
        if (i_next == pintervals->data.end())
        {
            DEBUG_CERR("    increment_flanking_interval_pair: Beyond any intervals after increment\n")
            set_next_out_of_bounds();
			return;
        }

		DEBUG_CERR("    increment_flanking_interval_pair: Between intervals after increment\n")
		return;
    }

    // 
    // Flanking intervals.
    // Make them point to the same interval
    // 
    i_prev = i_next;
    DEBUG_CERR("    increment_flanking_interval_pair: Within single interval after increment\n")
    return;
}



//!_______________________________________________________________________________________
//!
//! 	t_interval_search_cache::contains
//! 
//!     	tests either within a single interval or
//!     	returns the closer of the two intervals
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param test_pos 		
//!     \param dist 				distance to closest interval
//!     \param closest_interval 	which interval is closer
//!     
//!     \return bool 				If false, neither dist or closest_interval is set
//!_______________________________________________________________________________________
bool t_interval_search_cache::contains(long test_pos, long& dist, ci_interval& closest_interval) const
{
	assert(pintervals);

	// 
	// Contained within single interval
	// 
	if (single())
	{
		if (test_pos >= i_prev->beg && 
			test_pos <  i_prev->end  )
		{
			dist = 0;
			closest_interval = i_prev;
			return true;
		}
		return false;
	}


	//
	//  Not contained between intervals
	//  	return false
	// 
	if (test_pos < i_prev->end ||
		test_pos >=  i_next->beg  )
		return false;



	long dist1 = test_pos - i_prev->end + 1;
	long dist2 = i_next->beg - test_pos;

	if (next_out_of_bounds())
	{
		// 
		//	both out of bounds: empty data: return false
		// 
		if (prev_out_of_bounds())
			return false;

		//
		//	Next out of bounds: return Prev 
		// 
		dist = dist1;
		closest_interval = i_prev;
		return true;
	}

	//
	//	Prev out of bounds or next closer: return Next
	// 
	if (prev_out_of_bounds() || dist1 > dist2)
	{
		dist = dist2;
		closest_interval = i_next;
		return true;
	}

	//
	//	Prev closer: return Next
	// 
	dist = dist1;
	closest_interval = i_prev;
	return true;
}




//!_______________________________________________________________________________________
//!
//! 	t_interval_search_cache::contains
//! 
//!     	tests either within a single interval or
//!     	returns the closer of the two intervals
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param test_pos 		
//!     
//!     \return bool 
//!_______________________________________________________________________________________
bool t_interval_search_cache::contains(long test_pos) const
{
	assert(pintervals);

	// 
	// Contained within single interval
	// 
	if (single())
		return (test_pos >= i_prev->beg && 
				test_pos <  i_prev->end  );

	// 
	//	both out of bounds: return false
	// 
	if (next_out_of_bounds() && prev_out_of_bounds())
		return false;

	// 
	// Between intervals?
	// 
	return(test_pos >= i_prev->end && 
		   test_pos <  i_next->beg  );
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     t_intervals search functions
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//!_______________________________________________________________________________________
//!
//! 	contains
//! 
//! 		whether an interval contains a search position
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param test_pos The position to search
//!     \param cached_interval_pair A cache to speed up the search
//!     
//!     \return bool 
//!_______________________________________________________________________________________
#ifndef	NDEBUG
// debug will automatically print out new interval
struct t_print_new_interval
{
	t_interval_search_cache& cached_interval_pair;
	t_print_new_interval(t_interval_search_cache& cached_interval_pair_):
			cached_interval_pair(cached_interval_pair_){}
	~t_print_new_interval()
	{cerr << "    New interval: " << cached_interval_pair.str() << "\n";}
};
#endif

bool t_intervals::contains(long test_pos, t_interval_search_cache& cached_interval_pair) const
{
	DEBUG_CERR("    Check if contains " << test_pos << "?\n")

	if (data.size() == 0)
		return false;

	//
	//	cache invalid
	// 
	if (cached_interval_pair.pintervals != this)
	{	
		DEBUG_CERR("    Different cached interval set from previous test\n")
		//
		//	Get next interval 
		// 
		get_flanking_intervals(test_pos, cached_interval_pair);
		return cached_interval_pair.single();
	}

	// 
	// if within previous flanking intervals.
	// Just test if within a single interval
	if (cached_interval_pair.contains(test_pos))
	{
		DEBUG_CERR("    Same interval as previous test: " << cached_interval_pair.str() << "\n")
		return cached_interval_pair.single();
	}

	#ifndef	NDEBUG
	cerr << "    Different interval. Previous = " << cached_interval_pair.str() << "\n";
	t_print_new_interval print_new_interval(cached_interval_pair);
	#endif


	// 
	// Test next interval(s) first
	// 
	unsigned num_to_test = cached_interval_pair.single() ? 2 : 1;

	for (unsigned i = 0; i < num_to_test; ++i)
	{
		++cached_interval_pair;
        if (cached_interval_pair.contains(test_pos))
			return cached_interval_pair.single();
	}

	//
	//	Return to proper search 
	// 
	get_flanking_intervals(test_pos, cached_interval_pair);
	return cached_interval_pair.single();
}

//!_______________________________________________________________________________________
//!
//! 	get_flanking_intervals
//! 
//!     	returns the two intervals which bracket a position
//!     	If the test position falls within the same position, the two intervals will be
//!     	the same, and cached_interval_pair will return single() == true
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param test_pos 
//!_______________________________________________________________________________________
struct t_sort_by_end
{
	bool operator()(const t_interval& me, const t_interval& other)
	{
		if (me.end == other.end)
			return me.beg < other.beg;
		return me.end < other.end;
	}
};
void t_intervals::get_flanking_intervals(long test_pos, t_interval_search_cache& cached_interval_pair) const
{
	// points to me now
	cached_interval_pair.pintervals = this;

	// 
	//	check not past the end 
	// 
	cached_interval_pair.i_next = lower_bound(data, t_interval(test_pos, test_pos), t_sort_by_end());


	// If test_pos is before first interval: 
	// 
	//  	prev_interval => out_of_bounds_intervals[0]
	//  	next_interval => intervals[0]
	// 
	if (cached_interval_pair.i_next == data.begin())
	{
		if (test_pos < cached_interval_pair.i_next->beg)
		{
			DEBUG_CERR("    get_flanking_intervals: Before any intervals\n");
			cached_interval_pair.set_prev_out_of_bounds();
			return;
		}
	
		// 
		// entirely contained:
		// both flanking intervals point to the same thing
		// 
		else
		{
			DEBUG_CERR("    get_flanking_intervals: Contained within a single interval\n")
			cached_interval_pair.i_prev = cached_interval_pair.i_next;
			return ;
		}
	}


	// 
	// after last interval
	// 
	//  	next_interval => out_of_bounds_intervals[1]
	//  	prev_interval => intervals[-1]
	// 
	if (cached_interval_pair.i_next == data.end())
	{
		DEBUG_CERR("    get_flanking_intervals: Beyond any intervals\n")
		// Between last interval
		cached_interval_pair.i_prev = data.begin() + data.size() - 1;
		// ... and out of bounds
		cached_interval_pair.set_next_out_of_bounds();
		return;
	}

	// 
	// entirely contained:
	// both flanking intervals point to the same thing
	// 
	if (test_pos >= cached_interval_pair.i_next->beg) // && test_pos < cached_interval_pair.i_next->end)
	{
		DEBUG_CERR("    get_flanking_intervals: Contained within a single interval\n")
		cached_interval_pair.i_prev = cached_interval_pair.i_next;
		return ;
	}

	// 
	// before next
	// 
	assert(cached_interval_pair.i_next != data.begin());
	cached_interval_pair.i_prev = cached_interval_pair.i_next - 1;
	DEBUG_CERR("    get_flanking_intervals: Between two intervals\n")
	return;
}

//!_______________________________________________________________________________________
//!
//! 	any_closest
//! 
//!     	returns the two intervals which bracket a position
//!     	If the test position falls within the same position, the two intervals will be
//!     	the same, and cached_interval_pair will return single() == true
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param test_pos 
//!_______________________________________________________________________________________
bool t_intervals::any_closest(long test_pos, long& dist, ci_interval& closest_interval, t_interval_search_cache& cached_interval_pair) const
{
	if (data.size() == 0)
		return false;

	//
	//	cache invalid
	// 
	if (cached_interval_pair.pintervals != this)
	{	
		DEBUG_CERR("    Different cached interval set from previous test\n")
		//
		//	Get next interval 
		// 
		get_flanking_intervals(test_pos, cached_interval_pair);
        cached_interval_pair.contains(test_pos, dist, closest_interval);
		return true;
	}

	// try current flanking interval pair
    if (cached_interval_pair.contains(test_pos, dist, closest_interval))
	{
		DEBUG_CERR("    Same interval as previous test: " << cached_interval_pair.str() << "\n")
		return true;
	}
	#ifndef	NDEBUG
	cerr << "    Different interval. Previous = " << cached_interval_pair.str() << "\n";
	t_print_new_interval print_new_interval(cached_interval_pair);
	#endif

	// or next one
    ++cached_interval_pair;
    if (cached_interval_pair.contains(test_pos, dist, closest_interval))
		return true;

	// or next one
    ++cached_interval_pair;
    if (cached_interval_pair.contains(test_pos, dist, closest_interval))
		return true;

	//
	//	OK. do proper search
	// 
	get_flanking_intervals(test_pos, cached_interval_pair);
    cached_interval_pair.contains(test_pos, dist, closest_interval);
	return true;
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     t_intervals functions
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//!_______________________________________________________________________________________
//!
//! 	len
//! 
//! 		Assumes overlapping combined 
//!     
//!     \author lg (1/7/2011)
//!     
//! 	\return long cummulative size of all itervals. 
//!_______________________________________________________________________________________ 
long t_intervals::len() const
{
	long total_interval_size = 0U;
	BOOST_FOREACH(const t_interval& interval, data) 
		total_interval_size += interval.end - interval.beg;
	return total_interval_size;
}












//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     Defined in terms of non-self-modifying
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//!_______________________________________________________________________________________
//!
//!     bracket_by
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param by how much to expand each interval by
//!_______________________________________________________________________________________
t_intervals& t_intervals::bracket_by(long by, bool do_not_combine_afterwards)
{
	t_intervals ret = bracketed_by(by, do_not_combine_afterwards);
	swap(ret);
	return *this;
}


//!_______________________________________________________________________________________
//! 
//!		complement
//!			
//! 		Side effect: sorts and combines
//! 
//!		\author lg (24/01/2011)
//!     
//!     \param first 
//!     \param last 
//!_______________________________________________________________________________________
t_intervals& t_intervals::complement(long first, long last)
{
	t_intervals ret = complemented(first, last);
	swap(ret);
	return *this;
}

//!_______________________________________________________________________________________
//!     
//! 	touch
//! 
//! 		Only keep intervals which touch those in other
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!_______________________________________________________________________________________
t_intervals& t_intervals::touch(t_intervals& other )
{
	t_intervals ret = touched_by(other);
	swap(ret);
	return *this;
}

//!_______________________________________________________________________________________
//!     
//!		intersect_with
//!			
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!_______________________________________________________________________________________
t_intervals& t_intervals::intersect_with(t_intervals& other )
{
	t_intervals ret = intersected_with(other);
	swap(ret);
	return *this;
}

//!_______________________________________________________________________________________
//!     
//!		subtract
//!			
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals 
//!_______________________________________________________________________________________

t_intervals& t_intervals::subtract(t_intervals& other )
{
	t_intervals ret = subtracted(other);
	swap(ret);
	return *this;
}

//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     non-self-modifying
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
t_intervals t_intervals::bracketed_by(long by, bool do_not_combine_afterwards) const
{
	if (data.size() == 0)
		return t_intervals();

	t_intervals ret;
	ret.data.reserve(data.size() * 2);

	ci_interval src = data.begin();
	ci_interval end = data.end();

	while (src != end)
	{
		// make sure no under/overflow
		ret.data.push_back(t_interval(max(by, src->beg) - by, src->beg, src->id));
		ret.data.push_back(t_interval(	src->end, 
										min(numeric_limits<long>::max() - by, src->end) + by, 
										src->id));
		++src;
	}
	if (!do_not_combine_afterwards)
		ret.combine_overlapping();
	return ret;
}


//!_______________________________________________________________________________________
//! 
//!		complemented
//!			
//! 		Side effect: sorts and combines
//! 
//!		\author lg (24/01/2011)
//!     
//!     \param first 
//!     \param last 
//!_______________________________________________________________________________________
t_intervals t_intervals::complemented(long first, long last)
{
	//
	//	make sure no empty range
	//
	if (first >= last)
	{
		data.clear();
		return t_intervals();
	}

	//
	//	sort and combine
	//
	combine_overlapping ();

	t_vec_intervals	new_data;
	ci_interval	ii  = data.begin();
	ci_interval	end = data.end();

	// 
	// skip past beginning
	// 
	//  	Note: we could do a binary search here but pay off too small
	//
	ci_interval	prev = ii;
	while (ii != end)
	{
		if (ii->beg > first)
			break;
		prev = ii;
		++ii;
	}

	//
	//	adjust first if first within an interval
	//
	long beg = first;
	if (prev != ii && ii != end)
		beg = max(first, prev->end);



	while (ii != end && ii->beg < last)
	{
		new_data.push_back(t_interval(beg, ii->beg));
		beg = ii->end;
		++ii;
	}

	//
	//	add last intervals if necessary
	//
	if (beg < last)
		new_data.push_back(t_interval(beg, last));


	t_intervals ret;
	ret.data.swap(new_data);
	return ret;

// 
//	equivalent python code 
// 
//  from collections import namedtuple
//  t_interval = namedtuple("t_interval", "beg end")
//  data = [                                                     
//          t_interval(0,10 ),                                   
//          t_interval(20,30),                                   
//          t_interval(31,40),                                   
//          t_interval(50,60),                                   
//      ]                                                        
//                                                               
//  def func(data, first, last):                                 
//      ii = 0                                                   
//      prev = ii                                                
//      while (ii < len(data)):                                  
//          if (data[ii].beg > first):                           
//              break                                            
//          prev = ii                                            
//          ii+= 1                                               
//      if (ii == len(data)):                                    
//          new_data = [t_interval(first, last)]
//          print new_data
//          return
//      print prev, ii                                           
//      new_data = []                                            
//      #                                                        
//      #   adjust first if fast within an interval              
//      #                                                        
//      beg = first                                              
//      if (prev != ii):                                         
//          beg = max(first, data[prev].end)                     
//      while (ii < len(data) and data[ii].beg < last):          
//          new_data.append(t_interval(beg, data[ii].beg))       
//          beg = data[ii].end                                   
//          ii += 1                                              
//      #                                                        
//      #   add last intervals if necessary                      
//      #                                                        
//      if (beg < last):                                         
//          new_data.append(t_interval(beg, last))               
//      print data                                               
//      print new_data                                           
//  
//  func(data, 0, 100)
//  func(data, 100, 200)
//  func(data, -100, -50)
}




//!_______________________________________________________________________________________
//!     
//!		intersected_with
//!			
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!_______________________________________________________________________________________
t_intervals t_intervals::intersected_with(t_intervals& other )
{
	// 
	// empty ranges intersect to nothing
	//
	if (!data.size() || ! other.size())
	{
		return t_intervals();
	}
	combine_overlapping ();
	other.combine_overlapping ();

	t_vec_intervals		new_data;

	ci_interval	ii		= data.begin();
	ci_interval	ii_end	= data.end();
	ci_interval	jj		= other.data.begin();
	ci_interval	jj_end	= other.data.end();

	while (ii != ii_end && jj != jj_end)
	{
		//
		//	No overlap: keep going 
		// 
		if      (ii->end <= jj->beg)
			++ii;
		else if (jj->end <= ii->beg)
			++jj;
		else
		{
			new_data.push_back(t_interval(	max(ii->beg, jj->beg), 
											min(ii->end, jj->end), ii->id));
			if (ii->end < jj->end)
				++ii;
			else if (jj->end < ii->end)
				++jj;
			else
			{
				++ii;
				++jj;
			}
		}
	}

	t_intervals ret;
	ret.data.swap(new_data);
	return ret;
}

//!_______________________________________________________________________________________
//!     
//!		touched_by
//!     
//!		\author lg (17 June 2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals 
//!_______________________________________________________________________________________
t_intervals t_intervals::touched_by(t_intervals& other ) const
{
	// 
	// empty ranges intersect to nothing
	//
	if (!data.size() || ! other.size())
	{
		return t_intervals();
	}
	other.combine_overlapping ();

	t_vec_intervals		new_data;

	ci_interval	ii		= data.begin();
	ci_interval	ii_end	= data.end();
	ci_interval	jj		= other.data.begin();
	ci_interval	jj_end	= other.data.end();

	while (ii != ii_end && jj != jj_end)
	{
		//
		//	No overlap: keep going 
		// 
		if      (ii->end <= jj->beg)
			++ii;
		else if (jj->end <= ii->beg)
			++jj;
		else
		{
			new_data.push_back(t_interval(ii->beg, ii->end, ii->id));
			++ii;
		}
	}

	t_intervals ret;
	ret.data.swap(new_data);
	return ret;
}



//!_______________________________________________________________________________________
//!     
//!		subtracted
//!			
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals 
//!_______________________________________________________________________________________
t_intervals t_intervals::subtracted(t_intervals& other )
{
	if (data.size() == 0)
		return t_intervals();
	if (other.size() == 0)
		return t_intervals(*this);

	t_vec_intervals	new_data;

	combine_overlapping();
	other.combine_overlapping();

	ci_interval	ii		= data.begin();
	ci_interval	ii_end	= data.end();
	ci_interval	jj		= other.data.begin();
	ci_interval	jj_end	= other.data.end();


	// can do this because we know that ii != data.end()
	long beg = ii->beg;
	while (ii != ii_end && jj != jj_end)
	{
		//
		//  jj behind:
		//  	keep going 
		// 
		if 	(jj->end <= ii->beg)
		{
			++jj;
			continue;
		}

		// 
		//  ii behind:
		//  	save tail end of ii before incrementing
		// 
		if  (ii->end <= jj->beg)
		{
			if (ii->end != beg)
				new_data.push_back(t_interval(	beg, 
												ii->end,
												ii->id));
			++ii;
			beg = ii->beg;
			continue;
		}


		//
		//  overlap
		// 		save any remainder at beginning of ii
		// 
		if (beg < jj->beg)
			new_data.push_back(t_interval(	beg, 
											jj->beg,
											ii->id));

		// 
		// ii longer than jj
		//  	increment jj
		//  	and set next segment to save (beg) from end of jj
		// 
		if (ii-> end > jj->end)
		{
			beg = jj->end;
			++jj;
			continue;
		}

		// 
		// matching ends
		// 		increment both
		if (ii->end == jj->end)
			++jj;

		//	
		//	set beg after increment ii
		++ii;
		beg = ii->beg;
	}

	if (ii != ii_end)
	{
		if (ii->end != beg)
			new_data.push_back(t_interval(	beg, 
											ii->end,
											ii->id));
		++ii;
		while (ii != ii_end)
		{
			new_data.push_back(t_interval(	ii->beg, 
											ii->end,
											ii->id));
			++ii;
		}

	}


	t_intervals ret;
	ret.data.swap(new_data);
	return ret;
}










//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     Defined in terms of self-modifying
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//!_______________________________________________________________________________________
//!
//!     expanded_by
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param by how much to expand each interval by
//!_______________________________________________________________________________________
t_intervals t_intervals::expanded_by(long by, bool do_not_combine_afterwards) const
{
	t_intervals ret(*this);
	ret.expand_by(by, do_not_combine_afterwards);
	return ret;
}



//!_______________________________________________________________________________________
//!
//! 	combine_close_neighbours
//! 
//!     	 Merge neighbours less than min_distance distance apart.
//!     
//!     \author lg (1/7/2011)
//!     
//! 	\param min_distance
//!_______________________________________________________________________________________
t_intervals t_intervals::close_neighbours_combined(long min_distance) const
{
	t_intervals ret(*this);
	ret.combine_close_neighbours(min_distance);
	return ret;
}

//!_______________________________________________________________________________________
//!
//! 	removed_intervals_smaller_than
//! 
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param min_length intervals smaller than this will be removed
//!_______________________________________________________________________________________ 
t_intervals t_intervals::removed_intervals_smaller_than(long min_length) const
{
	t_intervals ret(*this);
	ret.remove_intervals_smaller_than(min_length);
	return ret;
}

//!_______________________________________________________________________________________
//!     
//!		combined_with
//!			
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals 
//!_______________________________________________________________________________________

t_intervals t_intervals::combined_with(const t_intervals& other) const
{
	t_intervals ret(*this);
	ret.combine_with(other);
	return ret;
}


//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     self-modifying
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//!_______________________________________________________________________________________
//!
//!     expand_by
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param by how much to expand each interval by
//!_______________________________________________________________________________________
t_intervals& t_intervals::expand_by(long by, bool do_not_combine_afterwards)
{
	i_interval src = data.begin();
	i_interval end = data.end();

	if (src == end)
		return *this;

	DEBUG_CERR("Expand Orig : " << src->str() << " by " <<  by << "\n")
	while (src != end)
	{
		// make sure no under/overflow
		src->beg = max(by, src->beg) - by;
		src->end = min(numeric_limits<long>::max() - by, src->end) + by;
		++src;
	}

	DEBUG_CERR("New Interval: " << src->str() << "\n")

	if (!do_not_combine_afterwards)
		combine_overlapping();
	return *this;
}

//!_______________________________________________________________________________________
//!     
//!		sort
//!     
//!		\author lg (5/10/2011)
//!_______________________________________________________________________________________
void t_intervals::sort()
{
	if (unsorted)
	{
		using boost::range::sort;
		sort(data);
		unsorted = false;
	}
}


//!_______________________________________________________________________________________
//!
//! 	combine_close_neighbours
//! 
//!     	 Merge neighbours less than min_distance distance apart.
//!     	 Side effect:
//!     					1) Sort intervals
//!     					2) Invalidate cache
//! 						3) Keep id of first interval in set of overlapping
//!     
//!     \author lg (1/7/2011)
//!     
//! 	\param min_distance
//!_______________________________________________________________________________________
t_intervals& t_intervals::combine_close_neighbours(long min_distance)
{
	// ignore if no intervals to combine
	if (size() < 2)
		return *this;


	//
	//  sort
	// 
	sort();

	i_interval dst = data.begin();
	// can do this because there are at least 2 items
	i_interval src = data.begin() + 1;
	i_interval end = data.end();

	//
	//  merge in place, 
	// 
	while (src != end)
	{
		if (src->beg  < min_distance + dst->end + 1)
		{
			// merge into current dst
			dst->end = max(src->end, dst->end);
		}
		else
		{
			// keep current dst, inc and copy src->dst
			++dst;
			*dst = *src;
		}
		++src;
	}
	assert(dst != end);
	++ dst;
	data.erase(dst, end);
	return *this;
}




//!_______________________________________________________________________________________
//!
//! 	remove_intervals_smaller_than
//! 
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param min_length intervals smaller than this will be removed
//!_______________________________________________________________________________________ 
struct t_test_interval_size
{
	long min_length;
	t_test_interval_size(long min_length_):min_length(min_length_){}
	bool operator()(const t_interval& i) const { return (i.end - i.beg) < min_length;}
};

t_intervals& t_intervals::remove_intervals_smaller_than(long min_length)
{
	data.erase(	std::remove_if(	data.begin(), data.end(), 
								t_test_interval_size(min_length)),
				data.end());
	return *this;
}


//!_______________________________________________________________________________________
//!
//! 	combined_with
//! 
//! 		insert all of other, sort and combine
//!     
//!     \author lg (1/7/2011)
//!     
//!     \param other 
//!_______________________________________________________________________________________
t_intervals& t_intervals::combine_with(const t_intervals& other)
{
	data.insert(data.end(), other.data.begin(), other.data.end());
	unsorted = true;
	combine_overlapping();
	return *this;
}











//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     t_intervals_per_contig definitions
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888 

//!_______________________________________________________________________________________
//!     
//!		sort
//!     
//!		\author lg (5/10/2011)
//!_______________________________________________________________________________________
void t_intervals_per_contig::sort()
{
	for (auto ii = data.begin(), end = data.end(); ii != end; ++ii)
		ii->second.sort();
}


//!_______________________________________________________________________________________
//!
//! 	bracket_by
//! 
//! 		returns the two sides of the intervals only
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param by 
//!     \param do_not_combine_afterwards 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::bracket_by(long by, bool do_not_combine_afterwards)
{
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		ii.second.bracket_by(by, do_not_combine_afterwards);
	return *this;
}

//!_______________________________________________________________________________________
//!
//!     combine_with
//!     
//!     \author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig& 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::combine_with(const t_intervals_per_contig& other)
{
	//
	//  combine contigs in common
	// 
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		if (other.data.count(ii.first))
			ii.second.combine_with(other.data.find(ii.first)->second);
	//
	//  add contigs unique to other
	// 
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, other.data) 
		if (!data.count(ii.first))
			data[ii.first] = ii.second;
	return *this;
}


//!_______________________________________________________________________________________
//!
//!     combine_close_neighbours
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param min_distance 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::combine_close_neighbours(long min_distance)
{
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		ii.second.combine_close_neighbours(min_distance);
	return *this;
}


//!_______________________________________________________________________________________
//!
//!     complement
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param by 
//!     \param do_not_combine_afterwards 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::complement(long first, long last)
{
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		ii.second.complement(first, last);
	return *this;
}

//!_______________________________________________________________________________________
//!
//!     expand_by
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param by 
//!     \param do_not_combine_afterwards 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::expand_by(long by, bool do_not_combine_afterwards)
{
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		ii.second.expand_by(by, do_not_combine_afterwards);
	return *this;
}
 


//!_______________________________________________________________________________________
//!     
//!		touch
//!
//! 		Only keep intervals which touch those in other
//!     
//!		\author lg (17 June 2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig& 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::touch(t_intervals_per_contig& other )
{
	// 
	// 	remove all contigs not in other
	// 
	i_map_str_intervals ii 	= data.begin();
	i_map_str_intervals end = data.end();
	while (ii != end)
		if (!other.data.count(ii->first))
			data.erase(ii++);
		else
			++ii;

	// 
	//	intersect remaining contigs in common
	// 
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
	{
		assert(other.data.count(ii.first));
		ii.second.touch(other.data[ii.first]);
	}
	return *this;
}

//!_______________________________________________________________________________________
//!     
//!		intersect_with
//!			
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig& 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::intersect_with(t_intervals_per_contig& other )
{
	// 
	// 	remove all contigs not in other
	// 
	i_map_str_intervals ii 	= data.begin();
	i_map_str_intervals end = data.end();
	while (ii != end)
		if (!other.data.count(ii->first))
			data.erase(ii++);
		else
			++ii;

	// 
	//	intersect remaining contigs in common
	// 
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
	{
		assert(other.data.count(ii.first));
		ii.second.intersect_with(other.data[ii.first]);
	}
	return *this;
}


//!_______________________________________________________________________________________
//!
//!     remove_intervals_smaller_than
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param min_length 
//!_______________________________________________________________________________________
t_intervals_per_contig& t_intervals_per_contig::remove_intervals_smaller_than(long min_length)
{
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		ii.second.remove_intervals_smaller_than(min_length);
	return *this;
}

//!_______________________________________________________________________________________
//!     
//!		subtract
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig& 
//!_______________________________________________________________________________________

t_intervals_per_contig& t_intervals_per_contig::subtract(t_intervals_per_contig& other )
{
	// only subtact if other has contact
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		if (other.data.count(ii.first))
			ii.second.subtract(other.data[ii.first]);
	return *this;
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     non-self-modifying
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888


//!_______________________________________________________________________________________
//!
//! 	bracketed_by
//! 
//! 		returns the two sides of the intervals only
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param by 
//!     \param do_not_combine_afterwards 
//!_______________________________________________________________________________________
t_intervals_per_contig t_intervals_per_contig::bracketed_by(long by, bool do_not_combine_afterwards) const
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		ret.data[ii.first] = ii.second.bracketed_by(by, do_not_combine_afterwards);
	return ret;
}


//!_______________________________________________________________________________________
//!     
//!		combined_with
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig 
//!_______________________________________________________________________________________

t_intervals_per_contig t_intervals_per_contig::combined_with(const t_intervals_per_contig& other) const
{
	t_intervals_per_contig ret;
	//
	//  combine contigs in common
	// 
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		if (other.data.count(ii.first))
			ret.data[ii.first] = ii.second.combined_with(other.data.find(ii.first)->second);
		else
			ret.data[ii.first] = ii.second;
	return ret;
	
	//
	//  add contigs unique to other
	// 
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, other.data) 
		if (!data.count(ii.first))
			ret.data[ii.first] = ii.second;
	return *this;
}
//!_______________________________________________________________________________________
//!
//!     close_neighbours_combined
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param min_distance 
//!_______________________________________________________________________________________
t_intervals_per_contig t_intervals_per_contig::close_neighbours_combined(long min_distance) const
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		ret.data[ii.first] = ii.second.close_neighbours_combined(min_distance);
	return ret;
}

//!_______________________________________________________________________________________
//!
//!     complemented
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param by 
//!     \param do_not_combine_afterwards 
//!_______________________________________________________________________________________
t_intervals_per_contig t_intervals_per_contig::complemented(long first, long last)
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		ret.data[ii.first] = ii.second.complemented(first, last);
	return ret;
}

//!_______________________________________________________________________________________
//!
//!     expanded_by
//!     
//!     \author lg (1/25/2011)
//!     
//!     \param by 
//!     \param do_not_combine_afterwards 
//!     
//!     \return t_intervals_per_contig 
//!_______________________________________________________________________________________
t_intervals_per_contig t_intervals_per_contig::expanded_by(long by, bool do_not_combine_afterwards) const
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		ret.data[ii.first] = ii.second.expanded_by(by, do_not_combine_afterwards);
	return ret;
}


//!_______________________________________________________________________________________
//!     
//!		intersected_with
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig 
//!_______________________________________________________________________________________
t_intervals_per_contig  	t_intervals_per_contig::intersected_with(t_intervals_per_contig& other )
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		if (other.data.count(ii.first))
			ret.data[ii.first] = ii.second.intersected_with(other.data[ii.first]);
	return ret;
}

//!_______________________________________________________________________________________
//!     
//!		touched_by
//!     
//!		\author lg (17 June 2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig 
//!_______________________________________________________________________________________
t_intervals_per_contig  	t_intervals_per_contig::touched_by(t_intervals_per_contig& other ) const
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		if (other.data.count(ii.first))
			ret.data[ii.first] = ii.second.touched_by(other.data[ii.first]);
	return ret;
}

//!_______________________________________________________________________________________
//!
//!     removed_intervals_smaller_than
//!     
//!     \author lg (1/25/2011)
//!     
//!     \param min_length 
//!     
//!     \return t_intervals_per_contig 
//!_______________________________________________________________________________________
t_intervals_per_contig  	t_intervals_per_contig::removed_intervals_smaller_than(long min_length) const
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		ret.data[ii.first] = ii.second.removed_intervals_smaller_than(min_length);
	return ret;
}

//!_______________________________________________________________________________________
//!     
//!		subtracted
//!     
//!		\author lg (1/25/2011)
//!     
//!     \param other 
//!     
//!     \return t_intervals_per_contig& 
//!_______________________________________________________________________________________

t_intervals_per_contig t_intervals_per_contig::subtracted(t_intervals_per_contig& other )
{
	t_intervals_per_contig ret;
	BOOST_FOREACH(t_map_str_intervals::value_type& ii, data) 
		if (other.data.count(ii.first))
			ret.data[ii.first] = ii.second.subtracted(other.data[ii.first]);
		else
			ret.data[ii.first] = ii.second;
	return ret;
}



//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//
//     misc and search functions
//     
//8888888888888888888888888888888888888888888888888888888888888888888888888888888888888888
//!_______________________________________________________________________________________
//!
//!     get_sorted_contigs
//!     
//!     \author lg (1/10/2011)
//!     
//!     \return std::vector<std::string> 
//!_______________________________________________________________________________________
vector<string> 	t_intervals_per_contig::get_sorted_contigs(eSORT_CONTIGS sort_mode) const
{
	using boost::range::sort;
	using boost::range::copy;
	vector<string> sorted_contigs;
	copy(data | map_keys, back_inserter(sorted_contigs));

	if (sort_mode == eSORT_BY_EMBEDDED_NUMBERS)
		sort(sorted_contigs, t_sort_by_embedded_numbers());
	else if (sort_mode == eSORT_BY_CANONICAL_CHROMOSOMES)
		sort(sorted_contigs, t_sort_by_canonical_chromosomes());
	else
		sort(sorted_contigs);
	return sorted_contigs;
}


//!_______________________________________________________________________________________
//!
//!     len
//!     
//!     \author lg (1/10/2011)
//!     
//!     \return long 
//!_______________________________________________________________________________________
long t_intervals_per_contig::len() const 
{
	long total = 0U;
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		total += ii.second.len();
	return total;
}

//!_______________________________________________________________________________________
//!
//!     size
//!     
//!     \author lg (1/10/2011)
//!     
//!     \return long 
//!_______________________________________________________________________________________
t_vec_intervals::size_type t_intervals_per_contig::size() const
{
	long total = 0U;
	BOOST_FOREACH(const t_map_str_intervals::value_type& ii, data) 
		total += ii.second.size();
	return total;
}

//!_______________________________________________________________________________________
//!     
//!		contains
//!			forwards the search to the right contig
//!     
//!		\author lg (1/11/2011)
//!     
//!     \param contig 
//!     \param test_pos 
//!     \param cache_for_search 
//!     
//!     \return bool 
//!_______________________________________________________________________________________

bool t_intervals_per_contig::contains(const std::string& contig, long test_pos, t_interval_search_cache& cache_for_search) const
{
	ci_map_str_intervals ii = data.find(contig);
	// no matching contig
	if (ii == data.end())
		return false;
	return ii->second.contains(test_pos, cache_for_search);
}


//!_______________________________________________________________________________________
//!
//!     any_closest
//!     
//!     \author lg (1/10/2011)
//!     
//!     \param test_pos 
//!     \param dist 
//!     \param closest_interval 
//!     \param cache_for_search 
//!     
//!     \return bool 
//!_______________________________________________________________________________________
bool t_intervals_per_contig::any_closest(const std::string& contig, long test_pos, long& dist, ci_interval& closest_interval, t_interval_search_cache& cache_for_search) const
{
	ci_map_str_intervals ii = data.find(contig);
	// no matching contig
	if (ii == data.end())
		return false;
	return ii->second.any_closest(test_pos, dist, closest_interval, cache_for_search);
}





#include <boost/spirit/include/phoenix_core.hpp>
namespace phoenix = boost::phoenix;
#include "qi_parse_error_handler.h"
phoenix::function<errorHandlerT> errorHandler_genomic_intervals = errorHandlerT("Genomic Intervals");




template <typename Iterator>
struct t_genomic_interval_grammar
: qi::grammar<Iterator, t_parsed_genomic_interval()>
{
    t_genomic_interval_grammar(eINTERVAL_FILE_FORMAT file_format, unsigned& curr_line_number_, 
								bool exception_on_error, string& error_msg)
    : t_genomic_interval_grammar::base_type(genomic_interval), curr_line_number(curr_line_number_)
    {
		// 
		//	Parsing grammar 
		// 

		contig		%= -(qi::lit("Chr") | qi::lit("chr")) > +qi::graph;
		beg			%= qi::int_;
		end			%= qi::int_;
		ignored		= qi::omit[+qi::graph];
		whitespace  = qi::space;

		if (file_format == eRAW_FORMAT)
			genomic_interval %=	 				contig	// chr  
								 > whitespace >	beg     // beg  
								 > whitespace >	end;    // end; 
		else // eGFF_FORMAT
			genomic_interval %=	 			contig	// chr  
								> '\t' >	ignored     // build
								> '\t' >	ignored     // type  
								> '\t' >	beg    		 // beg  
								> '\t' >	end;   		 // end; 

        // give names for debug
        contig      .name("Chromosome");
        beg         .name("begin position");
		end			.name("end position");
        ignored		.name("Ignored field");

		using namespace qi::labels;
        qi::on_error<qi::fail>                               
        (genomic_interval, errorHandler_genomic_intervals(_1, _2, _3, _4, 
				phoenix::ref(curr_line_number), phoenix::ref(error_msg),
				exception_on_error));
    }
	unsigned& curr_line_number;
	
    qi::rule<Iterator, std::string()>     			contig;
    qi::rule<Iterator, unsigned()>      			beg;
	qi::rule<Iterator, unsigned()>      			end;
	qi::rule<Iterator, void()>    					whitespace;
	qi::rule<Iterator, void()>    					ignored;
    qi::rule<Iterator, t_parsed_genomic_interval()>	genomic_interval;
};






//void t_discrete_intervals::init()
//{
//	//	
//	//  Sort intervals
//	// 
//	using boost::range::sort;
//	for (i_interval_by_contig ii = intervals_by_contig.begin(), end = intervals_by_contig.end(); ii != end; ++ii)
//		sort((*ii).second);
//}


//qi::rule<string::iterator> create_genomic_interval_grammar(eINTERVAL_FILE_FORMAT file_format)
//{
//	// 
//	//	Parsing grammar 
//	// 
//	using qi::graph;
//	using qi::space;
//	using qi::int_;
//	using qi::lit;
//	using qi::omit;
//	//  Chr1	11	T	?
//	typedef string::iterator ITER;
//	qi::rule<ITER> whitespace = qi::omit[qi::space];
//	if (file_format == eRAW_FORMAT)
//		return 	-(lit("Chr") | lit("chr")) > +graph			> whitespace > 			// chr
//				int_										> whitespace >   		// beg
//				int_;														   		// end;
//	else // eGFF_FORMAT
//		return 	-(lit("Chr") | lit("chr")) > +graph			> '\t' > 		// chr
//				omit[+graph]								> '\t' >   		// build
//				omit[+graph]								> '\t' >   		// type
//				int_		 								> '\t' >   		// beg
//				int_;												   		// end
//}

qi::rule<string::iterator> create_named_genomic_interval_grammar(eINTERVAL_FILE_FORMAT file_format)
{
	// 
	//	Parsing grammar 
	// 
	using qi::graph;
	using qi::int_;
	using qi::lit;
	using qi::omit;
	//  Chr1	11	T	?
	typedef string::iterator ITER;
	auto whitespace = qi::omit[qi::space];
	if (file_format == eRAW_FORMAT)
		return 	-(lit("Chr") | lit("chr")) > +graph			> whitespace > 			// chr
				int_										> whitespace >   		// beg
				int_										> whitespace >   		// end;
				+graph;																// name
	else // eGFF_FORMAT
		return 	-(lit("Chr") | lit("chr")) > +graph			> '\t' > 		// chr
				omit[+graph]								> '\t' >   		// build
				+graph										> '\t' >   		// type
				int_		 								> '\t' >   		// beg
				int_;												   		// end

}

//________________________________________________________________________________________
//
//  read_intervals
// 
//________________________________________________________________________________________
void read_intervals(std::istream& f_intervals, t_intervals_per_contig& intervals_by_contig, 
					bool zero_based_coordinates, eINTERVAL_FILE_FORMAT file_format, 
					bool throw_exception_on_error, std::ostream* error_strm)
{
	t_map_str_intervals& data = intervals_by_contig.data;

	// adjust values to convert from inclusive to open close coordinates
	long adjust_value = (zero_based_coordinates ? 0U : 1U);
	unsigned cnt_line = 0;

	//sp_genomic_interval_grammar genomic_interval_grammar = create_genomic_interval_grammar(cnt_line);
	string error_msg;
	t_genomic_interval_grammar<string::iterator> grammar(file_format, cnt_line, throw_exception_on_error, error_msg);

	//qi::rule<string::iterator> grammar = create_genomic_interval_grammar(file_format);




	string line;
	while (std::getline(f_intervals, line))
	{
		cnt_line ++;

		// ignore comments
		if (!line.length() || line[0] == '#')
			continue;

		//if (!parse_genomic_interval(genomic_interval_grammar, line, interval))
		//  break;
		t_parsed_genomic_interval parsed_genomic_interval;
		auto line_beg = line.begin();
		if (!qi::parse(line_beg, line.end(), grammar, parsed_genomic_interval) && error_strm)
			*error_strm << format("Warning: Unexpected line #%1$d\n%2$s\n") % cnt_line % error_msg;

		// convert to open / closed by adding 1
		data[parsed_genomic_interval.contig].data.push_back(t_interval(parsed_genomic_interval.beg_pos - adjust_value, parsed_genomic_interval.end_pos));
	}

	intervals_by_contig.sort();
	//using boost::range::sort;
	//for (i_map_str_intervals ii = data.begin(), end = data.end(); ii != end; ++ii)
	//{
	//	sort(ii->second.data);
	//	ii->second.unsorted = false;
	//}
}

//!_______________________________________________________________________________________
//!     
//!		read_named_intervals
//! 		Read intervals
//! 		The id of each interval is the offset into the vector of interval names
//!     
//!		\author lg (3/4/2011)
//!     
//!     \param f_intervals 
//!     \param intervals_by_contig 
//!     \param interval_names 
//!     \param zero_based_coordinates 
//!_______________________________________________________________________________________
void read_named_intervals(std::istream& f_intervals, t_intervals_per_contig& intervals_by_contig, 
							vector<string>& interval_names, bool zero_based_coordinates, eINTERVAL_FILE_FORMAT file_format)
{
	t_map_str_intervals& data = intervals_by_contig.data;

	// adjust values to convert from inclusive to open close coordinates
	long adjust_value = (zero_based_coordinates ? 0U : 1U);
	unsigned cnt_line = 0;

	//sp_genomic_interval_grammar genomic_interval_grammar = create_genomic_interval_grammar(cnt_line, true);
	qi::rule<string::iterator> grammar  = create_named_genomic_interval_grammar(file_format);



	unordered_map<string, long> map_interval_names_to_indices;

	string line;
	while (std::getline(f_intervals, line))
	{
		cnt_line ++;

		// ignore comments
		if (!line.length() || line[0] == '#')
			continue;

		string contig, name;
		unsigned beg, end;
		auto line_beg = line.begin();
		if (file_format == eRAW_FORMAT)
		{
			if (!qi::parse(line_beg, line.end(), grammar, contig, beg, end, name))
				throw std::logic_error((format("Unexpected line #%1$d [%2$s]") % cnt_line % line).str());
		}
		else
			if (!qi::parse(line_beg, line.end(), grammar, contig, name, beg, end))
				throw std::logic_error((format("Unexpected line #%1$d [%2$s]") % cnt_line % line).str());

		// save interval name
		// interval id points to the index of the string vector 
		long interval_index;
		if (!map_interval_names_to_indices.count(name))
		{
			interval_index = map_interval_names_to_indices.size();
			map_interval_names_to_indices[name] = interval_index;
			interval_names.push_back(name);
		}
		else
			interval_index = map_interval_names_to_indices[name];


		// convert to open / closed by adding 1
		data[contig].data.push_back(t_interval(beg - adjust_value, end, interval_index));
	}

	intervals_by_contig.sort();
	//using boost::range::sort;
	//for (i_map_str_intervals ii = data.begin(), end = data.end(); ii != end; ++ii)
	//{
	//	sort(ii->second.data);
	//	ii->second.unsorted = false;
	//}
}

//________________________________________________________________________________________
//
//  write_intervals
// 
// 		read into vector<string, t_discrete_intervals>
//________________________________________________________________________________________
void write_intervals(	ostream& 						f_intervals, 
						const t_intervals_per_contig& 	intervals, 
						const string&					interval_type_name,
						bool 							zero_based_coordinates,
						bool							ignore_type_name
						)
{
	// adjust values to convert from inclusive to open close coordinates
	long adjust_value = zero_based_coordinates ? 0U : 1U;

	// get sorted contig names
	vector<string> sorted_contigs = intervals.get_sorted_contigs();

	namespace karma = boost::spirit::karma;
	using boost::spirit::karma::uint_;
	using boost::spirit::karma::lit;
    using boost::spirit::karma::generate;

	std::ostreambuf_iterator<char> sink(f_intervals);

	if (ignore_type_name)
	{
		BOOST_FOREACH(const string& contig_name, sorted_contigs) 
			BOOST_FOREACH(const t_interval& interval, intervals.data.find(contig_name)->second.data) 
				generate(sink,
							(
								karma::string 
								<< lit('\t') << uint_
								<< lit('\t') << uint_
								<< lit('\n')
							),
							contig_name, interval.beg + adjust_value, interval.end
						);
	}
	else
	{
		BOOST_FOREACH(const string& contig_name, sorted_contigs) 
			BOOST_FOREACH(const t_interval& interval, intervals.data.find(contig_name)->second.data) 
				generate(sink,
							(
								karma::string 
								<< lit('\t')<< karma::string 
								<< lit('\t') << uint_
								<< lit('\t') << uint_
								<< lit('\n')
							),
							interval_type_name, contig_name, interval.beg + adjust_value, interval.end
						);
	}
}



/*
//________________________________________________________________________________________
//
//  write_intervals_by_type
// 
// 		read into vector<string, t_discrete_intervals>
//________________________________________________________________________________________
void write_intervals_by_type(	std::ostream& 					f_intervals, 
								std::map<string, long>& 	map_interval_type_to_index, 
								std::vector<string> 			interval_type_names, 
								t_vec_discrete_intervals& 		intervals_by_contig_by_type,
								bool 							zero_based_coordinates)
{
	using boost::range::sort;
	using boost::range::copy;

	sort(interval_type_names);
	BOOST_FOREACH(const string& interval_type_name, interval_type_names) 
	{
		long interval_id = map_interval_type_to_index.find(interval_type_name)->second;
		write_intervals(	f_intervals, 
							intervals_by_contig_by_type[interval_id], 
							zero_based_coordinates,
							interval_type_name + "\t");

	}
}


//________________________________________________________________________________________
//
//  read_intervals_by_type
// 
// 		read into vector<string, t_discrete_intervals>
//________________________________________________________________________________________
void read_intervals_by_type(std::istream& 						f_intervals, 
							t_vec_discrete_intervals& 			intervals_by_contig_by_type, 
							std::map<string, long>& 		map_interval_type_to_index, 
							std::vector<string>& 				additional_interval_types, 
							bool 								zero_based_coordinates)
{
	// adjust values to convert from inclusive to open close coordinates
	long adjust_value = zero_based_coordinates ? 0U : 1U;
	long cnt_line = 0;

	sp_genomic_interval_by_type_grammar genomic_interval_by_type_grammar = create_genomic_interval_by_type_grammar(cnt_line);

	// assign as many copies as necessary for each interval_type in map
	intervals_by_contig_by_type.resize(map_interval_type_to_index.size());

	string str;
	while (std::getline(f_intervals, str))
	{
		cnt_line ++;
		if (!str.length() || str[0] == '#')
			continue;

		t_genomic_interval_by_type interval;
		if (!parse_genomic_interval_by_type(genomic_interval_by_type_grammar, str, interval))
			break; 

		long interval_type_index = 0;
		// additional type
		if (!map_interval_type_to_index.count(interval.interval_type))
		{
			map_interval_type_to_index[interval.interval_type] = intervals_by_contig_by_type.size();
			intervals_by_contig_by_type.push_back(t_discrete_intervals());
			additional_interval_types.push_back(interval.interval_type);
		}
		interval_type_index = map_interval_type_to_index[interval.interval_type];

		// convert to open / closed by adding 1
		intervals_by_contig_by_type[interval_type_index].intervals_by_contig[interval.contig].push_back(t_interval(interval.beg - adjust_value, interval.end));
	}

	//	
	//  initialise all
	// 
	using boost::range::sort;
	for (i_vec_discrete_intervals ii = intervals_by_contig_by_type.begin(), end = intervals_by_contig_by_type.end(); ii != end; ++ii)
		ii->init();
}

*/

}

